Magnetic anomaly inversion through the novel barnacles mating optimization algorithm

Dealing with the ill-posed and non-unique nature of the non-linear geophysical inverse problem via local optimizers requires the use of some regularization methods, constraints, and prior information about the Earth's complex interior. Another difficulty is that the success of local search algorithms depends on a well-designed initial model located close to the parameter set providing the global minimum. On the other hand, global optimization and metaheuristic algorithms that have the ability to scan almost the entire model space do not need an assertive initial model. Thus, these approaches are increasingly incorporated into parameter estimation studies and are also gaining more popularity in the geophysical community. In this study we present the Barnacles Mating Optimizer (BMO), a recently proposed global optimizer motivated by the special mating behavior of barnacles, to interpret magnetic anomalies. This is the first example in the literature of BMO application to a geophysical inverse problem. After performing modal analyses and parameter tuning processes, BMO has been tested on simulated magnetic anomalies generated from hypothetical models and subsequently applied to three real anomalies that are chromite deposit, uranium deposit and Mesozoic dike. A second moving average (SMA) scheme to eliminate regional anomalies from observed anomalies has been examined and certified. Post-inversion uncertainty assessment analyses have been also implemented to understand the reliability of the solutions achieved. Moreover, BMO’s solutions for convergence rate, stability, robustness and accuracy have been compared with the solutions of the commonly used standard Particle Swarm Optimization (sPSO) algorithm. The results have shown that the BMO algorithm can scan the model parameter space more extensively without affecting its ability to consistently approach the unique global minimum in this presented inverse problem. We, therefore, recommend the use of competitive BMO in model parameter estimation studies performed with other geophysical methods.


Scientific Reports
| (2022) 12:22578 | https://doi.org/10.1038/s41598-022-26265-0 www.nature.com/scientificreports/ parameters of buried magnetized sources can be achieved through some challenging processes [33][34][35][36][37] . Local and global search algorithms can be applied for the inversion procedure 21,38 . Derivative-based local search algorithms using various regularization procedures have been routinely used for this task. Despite their very fast convergence characteristics, they have some serious disadvantages. These algorithms require prior information depending on geological conditions. Moreover, their success largely depends on the choice of initial guess and therefore they cannot scan the entire model space. They typically tend to reach a minimum in the vicinity of the initial model parameter set. Thus, these optimizers may reach any local minimum instead of the global minimum which is the deepest valley in the objective function surface topography. On the other hand, derivative-free global optimization and metaheuristic algorithms use a random walk search to reach the minimum in the given model parameter space. They can effectively scan the entire model space and therefore do not require a well-designed initial guess. Additionally, they have the ability of escaping local minima in the complex nature of the error function topography. Because of these advantages, the use of such nature-inspired algorithms in geophysical inverse problems has become very popular over the last decade. In fact, the first applications of these intelligence algorithms date back to the years before 1940, which is called the pre-theoretical period when the studies were informal 39 . The period between 1940 and 1980 is the early period, and these algorithms were formally studied during this period. In the method-centric period which is between 1980 and 2000 metaheuristic studies increased dramatically and numerous different approaches were introduced as an alternative tool to classical optimization algorithms 40 . The period from the 2000s to the present is called the framework-centric period, and the idea of describing computational intelligence algorithms as frameworks rather than methods has increased 39 . Standard Particle Swarm Optimization (sPSO) and Genetic Algorithm (GA) are commonly used in geophysical inversion 41 . Some comparative studies have verified that sPSO outperforms GA in terms of accurateness and better convergence characteristics for various problems [42][43][44][45] . However, it should be noted that there is no optimal metaheuristic algorithm for solving all types of inverse problems. Therefore, new problem-specific global optimization algorithms are still being developed. Accordingly, adaptations of new nature-inspired derivative-free algorithms for the inversion of magnetic anomalies have taken their place in the literature. In most cases, the outputs obtained were not compared with the results of another metaheuristic algorithm. Besides, in most studies, possible uncertainties in the model parameters obtained were not investigated. However, it is a well-known fact that global optimization and metaheuristic algorithms allow to perform uncertainty analysis which is an essential step in understanding the reliability of the solutions obtained. The difficulty of the inverse problem is increased when multiple-source structures are used to represent the resulting magnetic anomalies. Therefore, some of the anomaly peaks are not used in most studies, and the number of source structures is reduced and attempts are made to find a solution. In addition, it is an important deficiency that the possible regional magnetic anomaly effect is mostly ignored. Table 1 lists information about magnetic anomaly inversion studies carried out in recent years with global optimization algorithms and the presented study. It is clear that comparative tests, uncertainty analyses and regional effect investigations were not performed together in previous studies. In this study, we present a novel bio-inspired algorithm called Barnacles Mating Optimizer (BMO) for the inversion of magnetic anomalies. This optimization algorithm imitates the unique mating character of barnacles. Three theoretical cases with different scenarios were used to understand the proficiency of the proposed algorithm. Besides, three real data sets, including a chromite ore anomaly (India), an uranium ore anomaly (India) and an intrusive Mesozoic dyke anomaly (Brazil) were taken into consideration. Prior to the inversion experiments, we performed some modal analyses by mapping the surface topographies of the objective function for the model Table 1. Magnetic anomaly inversion studies carried out in recent years with global optimization algorithms and the presented study.

Algorithm
Inspiration Compared with Uncertainty analyses Regional analyses Magnetic anomaly www.nature.com/scientificreports/ parameter pairs. Thuswise the resolvabilities of model parameters and their effects on objective function were explored. Following the previous works 9, 46 we also investigated the regional background effect through the second moving average (SMA) technique. Consistency and credibility of the solutions obtained were assessed by performing some uncertainty appraisal analyses. Finally, the performances of the BMO and sPSO algorithms in our case were compared.

Methodology
Forward modeling. The magnetic anomaly equations for some idealized sources ( Fig. 1) such as a sphere 60,61 , an infinitely long horizontal cylinder 62 , a thin dyke 23,63,64 and a thin sheet 63 are given below, respectively.
In the definitions given above, T is the total field magnetic anomaly, x i is the observation point on a profile, z 0 (m) represents the depth to the center of the buried magnetic source (spheres and cylinders) and the top of dykes and thin sheets, θ denotes the effective magnetization angle (degree), q represents the shape factor (dimensionless), x 0 defines the origin of the anomaly or the horizontal center coordinate of the magnetic source (m), and K denotes the coefficient of amplitude (nT × m 2q−2 ), which is related to the model shape. The shape factors 2.5, 2, 1, and 1 are used for spheres, horizontal cylinders, thin dykes, and thin sheets, respectively 12 . On the basis of principle of superposition, composite magnetic anomalies of multiple-source are calculated easily using the following definition 12 .
(1)  www.nature.com/scientificreports/ where T j and P are the magnetic anomaly of j-th source and the number of sources, respectively.

Undesired regional effect and local interruption.
It is well-acknowledged that, due to the heterogeneity of the Earth's interior, the measured residual magnetic anomaly of a certain shallow source is mostly corrupted by undesired regional effects and some local interferences. The mathematical definitions of these effects are given as follows 9 : where T gen (x i ) denotes the composite anomaly, T re (x i ) is the undesired regional effect (c 0 -c 3 are predefined constants), and T local (x i ) represents the local perturbations originated from interfering sources.

SMA method.
A residual magnetic response of a shallow source can be calculated by removing the regional effect from the observed data set. A recent study 9 reported the efficiency of the SMA method to achieve this goal using the following definition: where R 2 (x) is the approximated residual magnetic anomaly, and s (m) denotes the window length. Typically, in order to provide a reasonable estimate of the true model parameters we can simply use the mean output produced by delimiting the SMA magnetic anomalies, which is controlled by s values (window length). 65 , which is used to simulate the special mating behaviors of barnacles for the optimization of some engineering problems. Barnacles are mostly hermaphrodites and live in shallow and tidal waters. These fantastic individuals can attach themselves temporarily to substratum or symbionts in the water. For this purpose, they mostly use whales, sea snakes or any other crustaceans. In order to survive, barnacle mating groups surround neighbors and potential competitors within reach of their penises. The reported 65 satisfactory results of BMO in 23 challenging benchmark functions and power system analyses motivated us to apply it to nonlinear geophysical problems. Applications with many different scenarios are presented to show how the proposed algorithm can reach the global optimum without suffering from the local optima entrapments and without being affected by the ill-posed nature of the magnetic inverse problem. The steps of the algorithm are given in brief as follows.

BMO. BMO is a novel algorithm
Initialization. In BMO, it is assumed that barnacles are the solutions 65 . Thus, in the inverse problem presented here barnacles can be treated as the model parameters (K, θ, x0, z0, q). The matrix of population X is defined as follows: where dim is the number of control variables (dim = 5 in this study) and N denotes the number of barnacles (population size). The X i j (the i-th control variable of the j-th barnacle. Here, i = 1, 2, 3... dim, j = 1, 2, 3... N) given above is subject to the upper and lower bounds [lb i , ub i ]. The appraisal of X is terminated primarily, and the sorting process based on the obtained fitness values is carried out to detect the best solution (K, θ, x 0 , z 0 , q) at the top of X.
Selection process. BMO uses a different system compared to other evolution-based systems such as Genetic algorithm 66 , Differential Evolution 67 etc. As the selection of two barnacles is relied on the length of their penises (pl) a simple case can be presented 65 to illustrate this special mating behavior assuming that the best solution is located at the top of X at a particular iteration and the maximum penis length of barnacles pl = 7. Therefore, barnacle #1 is only able to mate with one of barnacles #2 and #7 (selected barnacles are located within pl). Hence, the following simple mathematical forms to achieve the selection process is proposed 65 : where randperm(N) is a function that returns a row vector that contains a random permutation of the integers from 1 to N inclusive. The barnacle d and barnacle m are the parents to be mated and they should locate within pl. www.nature.com/scientificreports/ However, if one of them does not fulfill this requirement, the normal mating process is therefore postponed. The next generation is updated using the sperm casting process instead.
Reproduction. BMO's reproduction progression differs from other evolution-based algorithms. There is no explicit formula for developing the reproduction of barnacles, so BMO produces offspring based on the principle of the Hardy-Weinberg 68,69 . The following expression is given for this procedure: where p represents random numbers in the range of [0, 1], q = (1 − p), X i barnacle d and X i barnacle m are selected parents using Eq. (9). p and q denote the percentage of characteristics of the mating objects implanted in the next generation. The offspring, therefore, inherits the behaviors depending on the probability of the random number in the range of [0, 1].
In the BMO algorithm pl has an important impact on the exploitation (normal mating process) and exploration stages. The exploitation process occurs if the selected barnacle to be mated is within the penis length of the male barnacle. If not, the sperm cast is implemented and considered as an exploration process of BMO as given below: where rand() denotes random numbers between [0, 1]. Along with Eq. (11), the new offspring is formed by the female barnacle because the female barnacle retrieves the sperm from the water left by the other barnacles. As can be seen from the brief introduction to BMO, pl is the algorithm-based control parameter, which needs to be selected before employing BMO. The tuning process of pl will be discussed later. Figure 2 describes the basic workflow and pseudocode of BMO algorithm.
Objective function and stopping criterion. Basically, the BMO algorithm uses model parameters as search agents, and then quantitatively simulated specific matching behaviors guide these search agents to converge iteratively. Finally, if the root-mean-square-error (RMSE) is minimized to a predefined small value or the (10) www.nature.com/scientificreports/ iteration number reaches the predefined maximum number (Iter_max), the iteration cycle is terminated. Subsequently, the model parameter set (K, θ, x 0 , z 0 , q) obtained at the last iteration is considered to be the optimal solution. The definition of RMSE used in our case is given as follows: where T obs is the observed data, T cal denotes the data calculated from the model solution, x i is the i-th data measuring point along the profile aforementioned, and M represents the number of observed points. The relative error between the true and estimated parameter sets is determined by: where abs() is an absolute function, S true and S pre denote the true and the estimated model parameter set.

Synthetic applications
Model analysis of the objective function. The ill-posedness characteristics of the geophysical inverse problem makes the model parameter estimation process suspect and error-prone 70,71 . Hence, pre-and postinversion uncertainty appraisal studies play a crucial role in the detection of the ambiguities and consistencies in the obtained source parameters 36 . Commonly used pre-inversion analyses in geophysical inverse problems is the study of the modal type of the objective function used. These analyzes allow us to make a preliminary assessment of the mathematical solvability of the inverse problem of interest. We used the following definition as the objective function in the applications.
where T p=4 (x i ) is the simulated magnetic anomaly generated by a multiple-source model, where its true parameter set (K, θ, x 0 , z 0 , q) is presented in Table 2. The following definition were used to generate the synthetic data set: Synthetic magnetic anomaly was generated along a 400 m long profile with 10 m data intervals. We considered a range of values that are half and twice (50% perturbations) the true model parameter values to produce a wide search space for the modal analyses. Figure 3 shows the distributions of the cost function. The middle points of these surface topography maps locate the true values of the source parameter pairs (the global minimum). The distribution maps of θ − K, z 0 − K and z 0 − θ pairs clearly revealed that the global minima values are enclosed www.nature.com/scientificreports/ with almost rounded contours and no local ones are situated, which indicates the uni-modal feature 72 . In short, the parameter pairs are uncorrelated and can most likely be estimated individually by executing an efficient inversion code. However, some elongated valleys and basins with different flat bottoms representing the lowest error function regions displaying the multi-modal 73 characteristics were observed in other maps Fig. 3), which make the problem complicated, decrease the resolvability characteristics of the source parameters, and increase the uncertainties in the solutions. This means that parameter pairs are correlated with each other. It is well-known fact that banana-shaped contours in objective function surface topographies around the global minimum make precise estimations almost impossible. In our case, we did not observe this kind of topographic features. Thus, model parameters of the non-linear inverse problem presented here can most likely be resolved. Based on these findings, it is clear that the objective function yields the composite modality 74 since the uni-modality and multi-modality exist together regarding all model parameter pairs. Error surface topographies of the objective function used indicate the characteristics of non-linearity, high dimensionality and various shapes in paired model spaces, which means that estimations of the model parameters are challenging. However, these difficulties can be overcome with a powerful inversion algorithm which has sufficient ability to approach the global minimum as close as possible without compromising its robustness in order not to be captured by the local minima. www.nature.com/scientificreports/ Parameter tuning process. The efficiency of the global optimization algorithm largely depends on finetuned control parameters. However, the control parameters proposed by the developers for their own problems or reported ones in previous studies are commonly used for many types of optimization problems. It must be noted that there is no algorithm-based control parameter value that can be successful for all types of inverse problems. These optimal values may vary depending on the nature of inverse problems 49,70,75 . Hence, a vital step of applying a metaheuristic effectively is to appropriately determine the best user-defined control parameter value/s. As already mentioned, pl plays a key role in defining the exploitation (normal mating process) and exploration stages (sperm cast process) of BMO. To understand the effect of pl values on the solutions a composite noise-free magnetic anomaly due to a multiple-source model was inverted (model parameters are given in Table 2). The search space values used in modal analyses were considered. We used 30 independent runs to suppress the stochastic nature of the metaheuristics so that a more objective evaluation and estimation can be achieved. The optimization processes were performed using 140 iterations and N = 80. As the pl of the BMO varies successively (from 0 × N to 1 × N), the obtained mean RMSE values are displayed in Fig. 4a. Figure 4b demonstrates the variation of the standard deviations (stds) of the calculated RMSE values, indicating the uncertainty of inversion. The variations in the curves (Fig. 4a,b) obtained via various pl values clearly showed that an ill-defined pl value can reduce the performance of the BMO significantly in terms of affecting the solution correctness and the stableness of searching the true parameter values. Additionally, contrary to our expectations, both smaller (e.g., pl = 0 × N) and larger (e.g., pl = 1 × N) pl values caused insufficient performances. By this way, instead of making a random guess based on intuition, the importance of determining the optimum control parameter value/s for the problem of interest was revealed. By performing many trial-and-error tests, we determined that 0.65 is the optimum value for the pl.
Inversion studies. We used three synthetic cases with different scenarios to investigate the performance of BMO. The first case discusses a residual magnetic anomaly originated from a multiple-source model given in Table 2. In the second case, we added a regional background effect using a 3rd order polynomial and a local anomaly caused by a near-surface interference (spherical source) to the first case anomaly. In the last case, the second case anomaly was contaminated with a certain degree of random noise (40%). The solutions obtained through the BMO algorithm were compared with those of the sPSO algorithm that uses the suggested 9,51 control parameters. Besides, we investigated the performance of the SMA operator in reducing the effect of the regional background. All experiments were conducted on a Windows 10 operating system with Intel(R) Core (TM) i5-6300HQ CPU (2.30 GHz) and 3.8 GB of RAM.
Case 1. The synthetic magnetic anomaly generated and the causative sources are shown in Fig. 5. We performed 30 independent runs using 80 search agents and 140 iterations for both BMO and sPSO. Table 3   www.nature.com/scientificreports/ were compared with the residual magnetic anomaly (Fig. 6a). Figure 6b,c show the behaviors of the mean value and the std of the calculated RMSE values against the iteration number. Figure 6d shows the calculated relative errors of the model parameters recovered via BMO and sPSO approaches. As seen from Fig. 6 and Table 3, the application of BMO algorithm produced stable, accurate, and robust performance. Model parameters with larger error values were obtained with the sPSO application. On the other hand, sPSO showed a faster convergence rate when using large population of search agents. sPSO and BMO converged to global minimum within 20 and 60 iterations, respectively. Although sPSO showed a faster rate of convergence, the resulting RMSE values showed that BMO better achieved an adequate trade-off to balance the exploration and exploitation stages.
Case 2. In this experiment we added regional and local effects to the synthetically generated magnetic anomaly using the following expressions: where T 2 gen (x i ) denotes the composite anomaly, T re (x i ) is the regional background effect and T local (x i ) represents the local perturbation originated from a near-surface spherical interference. Figure 7 illustrates the subsurface magnetized sources and the resulting magnetic responses. BMO and sPSO were applied using 30 independent runs with 80 search agents and 140 iterations. The control parameters of both optimizers were kept unchanged, and applications were performed without using the SMA technique. Search space bounds consistent with the first case were used. After each iteration, the mean, std of RMSE values, and the relative errors of estimated parameters were recorded to observe the performances (Fig. 8). Table 3 stores the quantitative results of estimations. The run-time of sPSO and BMO algorithms were 69 s and 66 s, respectively. It is clear that the mean RMSE value increased significantly with the final iteration. Moreover, the error values between the calculated and the true model parameters of four causative sources increased too. It is clear that both algorithms underperformed due to the presence of external effects incorporated. Subsequently, we experienced the effectiveness of the SMA technique in reducing the regional effects. Figure 9 illustrates the calculated SMA magnetic anomalies (red Figure 5. Simulated magnetic anomaly of multiple-source and the causative subsurface model.   www.nature.com/scientificreports/ lines) for some s values (0.7, 1.3, 1.6, 1.9, 2.2, and 2.5 × dx, dx equals the data spacing interval). Blue and black lines represent the calculated responses obtained through BMO and sPSO, respectively and are generally in well agreement with the red lines except for the anomalies generated from the interfering shallow spheric source. In this scenario, the running times of the sPSO and BMO algorithms are equal to the product of the times given previously and the s-values used, namely 69 s and 66 s. Hence, we will solely present the computational time considering the circumstance of delineating the residual anomaly. Outputs of this synthetic experiment and the detailed inversion results are shown in Fig. 10 and Table 3, respectively. It is clear that the SMA technique proved useful in eliminating the regional effect from the composite magnetic anomaly and therefore both algorithms yielded relatively agreeable results. However, examining the results in detail, it is seen that the mathematical nature of simulating the special mating behaviors of barnacles yielded lower fitness errors, more accurate model parameter values, and higher inversion stability than sPSO.
Case 3. In order to make the work of both algorithms more difficult, a noise content of 40% was added to the magnetic anomaly of case 2 by using the definition below: where mean() is the average value of the input, rand 1 () and rand 2 () return an array containing pseudorandom values between [0, 1] of a given size. Using the same computation procedures, we inverted the synthetic magnetic anomaly through both algorithms. The SMA noisy anomaly and estimated SMA magnetic responses obtained by means of BMO and sPSO (blue and black solid lines) for different s values are demonstrated in Fig. 11. Figure 12 and Table 3 show the outputs of the optimization via both algorithm and detailed inversion results, respectively. Applications showed that the added 40% random noise only produced minor perturbations to the SMA residual anomalies (red lines in Fig. 11). Both algorithms produced satisfactory solutions, but the details revealed the superiority of BMO in terms of lower error values, more accurate model parameter values, and more stable optimization. www.nature.com/scientificreports/ Post-inversion uncertainty appraisal analysis of synthetic case. As mentioned previously in the modal analysis section, the optimization problem presented here is mostly unstable and error-prone because of the extreme complexity of escaping massive local minima while effectively exploiting the unique global minimum. Therefore, post-inversion uncertainty assessment analyses are crucial for understanding the credibility of the model parameters recovered. To practice the uncertainty appraisal studies, we used the solutions obtained in the first case experiment. Best model parameters obtained in each independent runs were listed and then sorted in ascending order of their misfit values. We assumed An as a variable in order to determine the number of sorted parameters. That we used to compute the final results by calculating the mean responses. We considered An as a variable to determine the number of parameters sorted. We then calculated the mean of each sorted parameter values. Lastly, the mean magnetic anomaly response was obtained with the mean model parameter set. Figure 13 exhibits that BMO showed low sensitivity to the increment of An, however, the performance of sPSO was significantly affected with the use of larger An values (pointed part A). Consistently, this phenomenon is correlated with the larger std values in the inversion results of sPSO. On the other hand, the superiority of the anomaly-fitting ability of BMO regardless the value of An is clearly observed in pointed part B (Fig. 13). These findings clearly revealed that BMO is more capable than sPSO in obtaining more robust solutions for this inverse problem presented.

Real data applications
The performance of BMO in the inversion of real magnetic anomalies was tested using three field cases from India and Brazil. These magnetic anomalies were considered as residual responses and were studied by some researchers using several data processing techniques. Based on this information, we initially treated these data www.nature.com/scientificreports/ sets as residual magnetic anomalies and applied the BMO algorithm using pl = 0.65 × N. We applied sPSO using the same algorithm-based control parameters to obtain outputs for the comparison tests. In the experiments, both algorithms were run independently 30 times using 80 search agents and 140 iterations. We then applied BMO procedure to SMA anomalies to understand if these anomalies are contaminated by regional background effect and/or local interference. Finally, the model parameter solutions obtained were interpreted in an integrated manner with the known geological data and/or previous geophysical estimates given in Table 4.
According to prior information, the ore bodies in this area are known to have a pod-type structure that is more reminiscent of a sphere-shaped source 12 . The magnetic anomaly response 12,76 of these magnetized sources is presented in Fig. 14a with pink circles. The data sampling interval is 40 m for this anomaly. This anomaly was studied previously by using Very Fast Simulated Annealing (VFSA) technique with various source assumptions such as single spherical-shaped body, cylinder, dyke, and sheet 12 . This magnetic anomaly was interpreted to originate from three spherical source bodies 12 . Model parameter solutions revealed in that study are given in Table 4. First, we considered this finding and applied the BMO and sPSO algorithms to estimate the model parameters.
Search space bounds for model parameters are given in Table 5. Figure 14a displays the comparison between the observed data and the reconstructed mean responses. Figure 14b shows the calculated std values of obtained responses at each station. The mean value and std of RMSE against the iterations for 30 independent runs are shown in Fig. 14c,d. Table 5 lists the obtained model solutions (K, θ, x 0 , z 0 , q) and RMSEs. The run-time of sPSO and BMO were 65 s and 59 s, respectively. The inversion results indicated that the assumed three spherical bodies approximate the sources effectively. Despite the faster convergence rate of sPSO, BMO outperformed in terms of response-fitting and stability. In the next step, we applied BMO to the SMA anomaly (s = 0.5, 0.75, and 1 × dx, dx = 40 m). Figure 17a displays the SMA results (red lines) for various s values and obtained responses (blue www.nature.com/scientificreports/ lines). Notably, model parameters obtained using the SMA technique (Table 5) are in well agreement with the results accomplished from considering the acquired data as a pure residual response. This finding indicates that the observed magnetic response does not contain undesired structural information such as regional background and near-surface source effects, but may contain random noise. The solutions obtained here are relatively in well agreement with the ones achieved via VFSA (Table 5).
Uranium Ore Anomaly, India. The magnetic anomaly sampled at 8 m interval in the Beldih mine (Purulia, West Bengal, India) is caused from uranium ore bodies 77 . Some previous studies revealed that the uranium ore bodies in this part can be represented with vertical thick-sheet-like structures as deciphered from spontaneouspotential data interpretation [78][79][80] . This interpretation was validated by drilling results 81 , which disclosed that mineralization begins from the near-surface and extends to the depths of 10-20 m and which are approximately vertical and dipping northerly to southerly 12 . Additionally, gravity and resistivity responses were also used to understand subsurface nature of the source structures 77,82,83 . Lastly, this magnetic anomaly was interpreted 12 using the VFSA technique considering dyke-like source as equal to thick-sheet-type source (Table 4). Similarly, we considered it as a residual response due to multiple dyke-like structures and implemented BMO and sPSO algorithms for reinterpretation. The model parameter ranges of assumed dykes and the results obtained are listed in Table 5. The run-time of sPSO and BMO were 54 s and 48 s, respectively. The comparison of the calculated magnetic responses with the observed one is illustrated in Fig. 15a. Figure 15b depicts the calculated std values  Figure 15c,d demonstrates the mean value and std values of RMSE against the iteration number. It can be clearly seen that the anomaly of the assumed multiple dyke-like structures matches well with the observed one and BMO provides appealing results in terms of escaping local minima, fitting errors and stability of the inversion compared to those obtained via sPSO. Figure 17b illustrates the SMA anomalies (red lines) for s = 0.5, 0.75, and 1 × dx (dx = 8 m) and reconstructed responses via BMO (blue lines). An obvious difference between the calculated and the reconstructed responses is seen. Nevertheless, recovered model parameters given in Table 5 are correspond to the solutions obtained without the implementation of the SMA scheme. Thus, it was verified that the observed anomaly is most likely corrupted by the local interference effect and/or random noise.
Mesozoic dyke Anomaly, Brazil. Figure 16a shows a magnetic anomaly observed on a Mesozoic dyke that intruded into Paleozoic sedimentary rocks in the Parnaiba basin, Brazil 84 . The data sampling interval is 0.8 m. This anomaly was previously interpreted using different approaches 12,26,84,85 assuming a horizontal cylinder along with thin-sheet structure. Solutions reported in those studies are given in Table 4. Here, this magnetic anomaly was studied using BMO and sPSO algorithms, considering it the response of a cylindrical body. Model parameter search bounds used are listed in Table 5. The run-time of sPSO and BMO were 20 s and 18 s, respectively. Estimated mean responses through BMO and sPSO are displayed in Fig. 16a, which closely match with the measured one (pink circles). However, the unsatisfactory performances of sPSO and BMO are observed at the pointed part shown in Fig. 16a. The possible explanation for this case is the presence of the local interference effect. Figure 16b demonstrates the calculated std values of obtained responses at each measurement station. Figure 16c,d display the mean value and std values of RMSE against the iterations. Table 5 gives the estimated model parameters. Findings showed that BMO yields better performance in terms of relatively lower misfit values and higher stability again. Figure 17c shows a comparison between the calculated SMA results for various s values (s = 0.5, 0.75, and 1 × dx, dx = 0.8 m), and regenerated responses using BMO. The estimated source parameters are listed in Table 4 along with misfits. Accordingly, it can be mentioned that BMO produced steady solutions www.nature.com/scientificreports/ with those delineated from considering the observed data as a pure residual response. Calculated anomalies of the BMO, however, revealed obvious differences most likely due to the existence of local interference effects and/or random noises. Thus, we concluded that the observed magnetic anomaly does not contain remarkable regional component, but may contain some amount of local perturbation and/or random noise.

Post-inversion uncertainty appraisal analyses of real data cases.
The reliability of the model parameters obtained was investigated, as in the synthetic experiments on the basis of the solutions of 30 independent runs for all real data cases. Figure 18 shows the variation characteristics of computed responses using different An values. It can be easily observed that the performance of the sPSO was affected with a larger An regarding the first two multiple-source model, especially in the pointed parts. This fact correlates with the large std values as discussed before. On the contrary, BMO is generally insensitive to a larger An, which indicates the superiority of the BMO in handling more complicated cases. Notably, as for the single structure model, only minor variations occurred in the responses of BMO and sPSO. However, when comparing the best performances of both approaches using An = 2 it is again apparent that BMO provided better fitting performance in these real data cases, since it offers sufficient opportunities to escape massive local minima while effectively exploiting the global minimum.

Conclusions
BMO, a novel nature-inspired derivative-free global optimizer, was presented for the inversion of magnetic anomalies. Modal analyses helped us understand the characteristic of the inverse problem. Finely-tuned control parameter increased the effectiveness of the optimizer. Satisfactory solutions were obtained in both synthetic and real data cases. SMA technique successfully suppressed the regional background effects. The effectiveness of BMO was compared with sPSO, the most widely used metaheuristic in geophysical inverse problems. The reliability of the results produced by both algorithms was tested with uncertainty analysis. Applications showed that BMO produced more stable and more robust solutions for the cases used in this study. It was also observed that BMO completed the optimizations in relatively less run-time than sPSO. Besides, the model parameters  www.nature.com/scientificreports/ obtained for the real data cases are in agreement with the available geologic and/or geophysical information. Based on these findings, we concluded that BMO is more efficient than sPSO in balancing global exploration (the ability to explore the search space extensively and escape the capture of massive local minima) and the local exploitation (intensively approximating the global minimum) processes in this type of inverse magnetic problem. The proposed optimizer quantitatively simulates the specific mating behavior of barnacles leads to a suitable trade-off between these two significant key processes. This study therefore represents an entirely new and competing tool for the efficient interpretation of causative sources from a variety of geophysical anomalies. On the other hand, similar to other metaheuristic algorithms, one of the weaknesses of the BMO is that it uses an algorithm-based control parameter (pl) whose optimum value may change according to each inverse problem. Therefore, parameter tuning studies are needed for each inverse problem in order to get the most out of the algorithm. Another weak point is that it requires more run-time than derivative-based local search algorithms. In addition, the increase in the number of model parameters to be estimated or the studies in which dense anomaly equations are used increase the run-times even more. However, the use of powerful and capable computers produced in line with technological developments can easily reduce this disadvantage.

Data availability
The data are available upon request from the first author.